Minimizing the population extinction risk by migration 
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Many populations in nature are fragmented: tliey consist of local populations occupying separate 
patches. A local population is prone to extinction due to the shot noise of birth and death processes. 
A migrating population from another patch can dramatically delay the extinction. What is the 
optimal migration rate that minimizes the extinction risk of the whole population? Here we answer 
this question for a connected network of model habitat patches with different carrying capacities. 
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Many populations in nature are fragmented. Such 
meta-populations consist of local populations occupying 
separate habitat patches [iHl]- Habitat fragmentation 
is implicated in the decline and extinction of many en- 
dangered species [1]. To mitigate the negative impact 
of habitat fragmentation, conservation biologists have 
called for the construction of corridors to facilitate mi- 
gration between separate habitat patches Q . Predicting 
how migration affects population persistence is impor- 
tant for species conservation, especially when the local 
population size is depressed, and the local populations 
become prone to extinction because of randomness of the 
birth and death processes. In this situation, it is of cru- 
cial importance to determine the optimal migration rate 
that maximizes the mean time to extinction (MTE) of 
the meta-population. This problem has attracted much 
of attention from ecologists, and has been addressed, for 
different meta-populations, in experiments and stochas- 
tic simulations [6l4ll . Here we approach this impor- 



tant problem theoretically, for a simple logistic model 
of stochastic local populations coupled by migration. We 
analyze rare large fluctuations causing population extinc- 
tion and show that there is an optimal migration rate 
that maximizes the MTE of the meta-population. 

Meta-population model. Mathematical biolo- 
gists have proposed different types of stochastic meta- 
population models. In a widely used class of models 
the local population distribution, its dynamics within a 
patch, and its effect on migration are ignored 
We show here that it is a proper account of these features 
that leads to the qualitatively new effect of the existence 
of an optimal migration rate. 

Consider N local populations of particles A located on 
a connected network of patches i = 1,2,...,A^. The par- 
ticles undergo branching A — >■ 2A with rate constant 1 
on each patch and annihilation 2A — >■ with rate con- 
stant l/{KiK) on patch i. The parameters Ki = 0{\), 
i — 1,2, . . . , N , describe the disparity among the local 
carrying capacities KiK . Each particle can also migrate 
between connected patches i and j with rate constant 
fiij — II ji. We assume that /iy = fJ.Mij, where elements 
of Mij are of order unity. 



For K ^ 1 each local population is expected to be 
long-lived. Still, the shot noise will ultimately drive the 
whole meta-population to extinction. The MTE of the 
meta-population, T, is exponentially large in K but fi- 
nite [nl, [l^l • How does T depend on the characteristic 
migration rate /x? At /i = each local population goes 
extinct separately, and 7),=o is determined by the patch 
with the greatest carrying capacity. Km = K inaxi{Ki}: 



lnTf,=o/K ~ 2(1 - In 2) max{Kj 



(1) 



(T for a single patch was found in Refs. [iMi)- The 
ideas behind our results for /i > are the following. At 
very fast migration, /i — >■ cx), the local populations be- 
come fully synchronized: both at the level of the expected 
local carrying capacities, and at the level of large fluctua- 
tions leading to population extinction. The total carrying 
capacity of the meta-population, as derived from the rate 
equation for this model (l^, becomes kK, where 



(2) 



One can argue, therefore, that at fj, ^ oo the meta- 
population goes extinct as if it were occupying a single 
effective patch with the total rescaled carrying capacity 
K, that is 
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lnr„/i^ ^2(l-ln2)K. 



(3) 



The main result of our work is that, for unequal Ki, T 
reaches its maximum at a finite value of the migration 
rate. This fact is intimately related to synchronization 
of the most probable local extinction events that occurs 
already at very small migration rates. The synchroniza- 
tion makes T close to that for a single patch with the 
combined carrying capacity k^: 



lnT^^o/A'^2(l-ln2)EK, 



(4) 



Now let us inspect the MTE as described by Eqs. ([T]), 
and As Ki > max^jKi} and Ki > k for any Ki, 
the MTE must reach a maximum at a finite value fj, = ii„, 
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unless all the patches have the same carrying capacity. 
We will present evidence that /x* <C 1 and scales as 1/K. 

How to understand qualitatively the non-trivial de- 
pendence of the MTE on jil Consider first the large- /i 
regime. Equation ^ implies that patches with smaller 
carrying capacities dominate the effective annihilation 
rate. For example, in a system of two patches, each par- 
ticle spends half its time on each of the two patches. 
Then the patch with the smaller carrying capacity dom- 
inates the total annihilation rate. As /i decreases, par- 
ticles will spend enough time on the good patch so that 
the total carrying capacity will drift up, and the MTE 
will increase. Now consider a very small but finite /i, so 
that the migration rate is higher than the (exponentially 
small) local extinction rates. Here, for the whole meta- 
population to go extinct, all local extinction events must 
occur in synchrony, and this leads to Eq. Q. 

Now we expose our results in more detail. For sim- 
plicity, we will first consider a system of two patches and 
then generalize our results to a network of N patches. 
The rate equations for the two-patch system are: 



y 



y 



K 



+ fJ'X - ny, 



(5) 



where x and y are the local population sizes rescaled 
by KiK, and k = K2/K1. Equations ^ have two 
fixed points: the unstable point Xq = yo = that de- 
scribes an empty system, and a stable point [a;*(K,/x) > 
0, y^,{K,fi) > 0] that describes an established meta- 
population. At ^ — one has x^ = 1 and = n, 
whereas for infinitely fast migration, /i — > 00, 



x^, ~ y^, ~ 2k/ (1 + k). 



(6) 



The characteristic time tr of population establishment 
is determined by the smaller of the two eigenvalues of 
the linear stability matrix of Eqs. ([5]) at the fixed point 
(x* , y* ) . 

In a stochastic formulation, the probability Pm,n{t) to 
find m particles in patch 1 and n particles in patch 2 
evolves in time according to the master equation 

Pm,n{t) = HPm.n = {m - l)Pm-l,n + [n - l)Pm,„_l 



(m + l)(m + 2) 
2K 



7n+2,n 



(n + l)(n + 2) 
2nK 



P 



7n,n+2 



+fi{in + l)P„i+i,,i-i + + l)P„i-i,„+i 
m(TO — 1) n{n — 1) 



(1 + fi){'m -I- n) -I- 



2K 



2kK 



Prr 



(7) 



The probability Pq.o that the meta-population goes ex- 
tinct by time t is governed by the equation 



PoAt) 



-P2,0 H T^Poa- 



(8) 



K ' kK 

Long-time dynamics and the MTE. For t > tr, 
Pm,n{t) becomes sharply peaked at the local carrying ca- 
pacities = if X* and = Ky.^, , corresponding to the 



stable fixed point (a;*,y*) of the mean- field theory. The 
subsequent slow decay of Pm.n in time is determined by 
the lowest excited eigenmode T:rn,n of the master equa- 
tion operator H: Pm,n{i) — 'I'm.n exp(— i/T). Simulta- 
neously, a probability peak at m = n = grows with 
time: Po,o(i) - l-exp(-t/T) The inverse eigen- 

value T is an accurate approximation to the MTE. Since 
it turns out to be exponentially large with respect to 
A' 3> 1, one can neglect the right-hand-side of the eigen- 
value problem Hnm.n = T^m.n/T and consider the quasi- 
stationary equation HTTm,n — 0. Once TTm.n is found, the 
MTE can be determined from Eq. ([5]): 



(9) 



WKB THEORY. To find iTm.n for not too small values of 
/i, we employ a dissipative variant of Wentzel-Kramers- 
Brillouin (WKB) approximation, pioneered in Refs. [25I- 
[28| . and extensively used in the problems of stochastic 



population extinction [17|, [181, [231, [24 
lii. The WKB ansatz is 



29-381, see also Ref. 



exp[-KS{x,y)], 



(10) 



where x = m/K and y = n/K are treated as continu- 
ous variables. We plug Eq. (|10p into the quasi-stationary 
equation H-Km.n — and Taylor expand S around {x,y). 
In leading order in 1/K <i; 1 this gives a zero-energy 
Hamilton- Jacobi equation H(x,y,dxS,dyS) — with 
classical Hamiltonian 

H{x.y.p„p,)=x{e<'- - 1) + y (e"*- - l) 



+y {ePy - 

+HX (e-p^+Py -i)+tiy {eP'^-Py - 1) . (11) 



The established population corresponds to the fixed point 
M — {x^,,y^.,0,0) of the Hamiltonian flow. Up to a 
pre-exponent, T ~ exp(KS), where S is the action 
along the instanton: a special zero-energy (H = 0) tra- 
jectory in the phase space {x,y,Pa:,Py) that exits, at 
time t — — oo, the flxed point M and approaches the 
fluctuational extinction point F that, for the two-patch 
branching-annihilation model, is (0, 0, — oo, — oo) [40|. In 
the absence of an independent integral of motion in ad- 
dition to the Hamiltonian itself, this trajectory, and the 
action along it, can only be found numerically. Analyti- 
cal results are possible in the limits of small and large /i 
that we will now consider. 

When fj, ^ the Hamiltonian (ITU) becomes separable, 
and the instanton trajectory can be easily found: 



x{t) = q{t - T^) , y{t) ^ Kq{t - Ty) 

Px{t) ^ p{t ~ T^), py{t)^p{t-Ty), (12) 



where 



g(t) = 2(2 + 3e'-f e^*)-\ p(t) = - ln(l + e*). (13) 
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Notice that the solution for /i — >■ includes arbitrary time 
shifts Tx and Ty in the x- and y populations, respectively. 
These will become important shortly. The action 

/oo 
[pxX +pyy - H) dt 
-oo 

= 2{l-\n2){l + K)c^lnT^^o/K. (14) 



Equation ([141) coincides with that for an effective one- 
patch system with the combined carrying capacity (1 + 
k)K. This extinction time is exponentially large com- 
pared with the one obtained if one neglects migration 
completely, see Eq. ([T]) with max{Ki} = 1. The sharp 
increase of T once slow migration is allowed results from 
synchronization of the most probable local extinction 
paths For /z <C 1, the two noisy local populations 

behave almost independently for typical, small fluctua- 
tions. For rare large fluctuations, such as the one causing 
extinction of the whole meta-population, the dynamics of 
the local populations becomes synchronized. How does 
the synchronization show up in the WKB calculations? 
In the absence of migration, fi = 0, the time shifts 
and Ty which appear in Eqs. (fT2)) are arbitrary, reflecting 
the time-translational invariance of local extinctions. A 
small /i > partially breaks this invariance and selects a 
particular relative time shift t ^ Ty ~ t^, implying syn- 
chronization. Since the zero-order action (jl4p is invari- 
ant with respect to the local time shifts, it is necessary 
to consider a small fj, correction, S — ^ 0) -f- AS 
in order to determine r. The first order correction can 
be calculated by integrating over the unperturbed x- and 
y-instantons (|12p : 



AS 



maxx(T), 



X{r) = / {q{t) 
+ Kq{t - t) 



-P{t)+p{t-r) 



- 1 



,p(t)-p(t-r) 



1 



I dt. (15) 



That is, the optimal time shift r = r, (k) is determined 
from the minimization of the action, that is the maxi- 
mization of x(''")j with respect to r. This minimization 
can be easily performed, as the integral in Eq. (jisp can 
be evaluated analytically (20| . 

By virtue of Eq. ([12]), x(0) = 0. This implies that 
AS < 0, and so T is a non-increasing function of ^ for 
/i <C 1. The function xi''') depicted in Fig. [T] for 
K — 1, 0.5 and 0.25. For k = 1 (two identical patches) 
the maximum is achieved at r = 0, as expected from 
symmetry, so AS* = 0. In this case the solution (fT2|) 
with K — 1 holds for all fi. That is, a higher migration 
rate docs not affect T up to a pre-exponential factor. 
For K < 1 we obtain Th.(k) > and AS < 0, that is T 
goes down with an increase of /i, see Fig. [2] Because 
of the large factor K, a small decrease in S translates 
into an exponentially large reduction of T of the meta- 
population. Note that the WKB approximation, leading 



to Eq. ([T5|) . is only valid for n ^ K~^. We expect that, 
for n < ^/K (but not exponentially small in K), weak 
synchronization [to within time uncertainty of {iJ,K)~^] 
occurs, again leading to MTE as in Eq. (|T4|) [i^ . 




FIG. 1; (Color online) Function x(''")) see Eq. (fTSt . for two 
patches, for k = 1 (dotted line), 0.5 (solid line) and 0.25 
(dashed line). 

Now consider the opposite limit, — >• oo. Here the 
total population size Q = x + y varies slowly in compar- 
ison with the fast migration. The fast variables x and y 
rapidly adjust to the slow dynamics of Q, staying close 
to their stationary values for the instantaneous value of 
Q. Transforming to Q and q ^ x and associated con- 
jugate momenta as a new set of canonical variables, one 
arrives [20] at a Hamiltonian, associated with the large 
fiuctuations of the total population size: 



H,ioUQ,P)^eHi[q{Q),Q,0,P 

1 + K 



1 + 



8k 



2 ^e-2P _ ^^ 



(16) 



with e = l//i ^ 1. Equation ([T5)) describes an effective 
single-patch Hamiltonian with a rescaled carrying capac- 
ity K = 4k/(1 + k), and we obtain 



lnT„ 



8(1 -ln2)K 



K 



1 



(17) 



For N — 2 this agrees with the announced result ([3]). 

WKB NUMERICS. For intermediate values of /z the 
instantons, and the associated action, can be found nu- 
merically: either by a shooting method [2^ H^], or by 
iterations 17,|33, |4l|. Here we used both methods, and 



the results for InT/K agreed within less than 1 per cent. 
Figure [2] shows the numerically found S for k = 0.25 and 
different /i, respectively. At fi <^ 1 the numerical results 
agree with the prediction of linear theory, Eq. (fT5|) . At 
large fj, they approach the asymptote ()17|) . Similar re- 
sults were obtained for other values of k. Figure [5] also 
compares the WKB results with those of a numerical so- 
lution of (a truncated version of) the full master equation 

Beyond WKB theory. To evaluate the maximum 
MTE and the optimal migration rate, one needs to re- 
solve the jump of {\nT)/K at /z = predicted by the 
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FIG. 2: (Color online) InT vs. the migration rate fi for 
two patches, n = 0.25. Circles: numerical WKB solutions. 
Diamond and square: predictions of Eqs. ([!]) and (|14p . re- 
spectively. Dashed line: prediction of Eq. (|15|) for ^ ^ 1: 
5S = -0.5455 Dotted line: prediction of Eq. dTT)) for 

2> 1. The solid line was obtained from a numerical solution 
of the master equation for K = 220. 



WKB theory, see Eqs. ([T|) and ([H]). We determined the 
MTE for exponentially small /i by numerically solving 
the master equation ([7]) and by performing stochastic 
simulations. The resulting /i-dependence of the MTE, at 
K = 0.25 and different K, is shown in Fig. |31 The maxi- 
mum of T is observed at a small migration rate fi^, that 
apparently scales as K~^. 





FIG. 3: (Color online) R-HnT vs. ^ (a) and In/i (b) for a 
two-patches meta-population from a numerical solution of the 
master equation and stochastic simulations, (a) k — 0.25 and 
K = 20,30,40 and 50 (bottom to top). Inset: the migration 
rate fit, at which the maximum of MTE is observed, versus 
K. (b) K = 0.25 and K = 20; dashed line: Eq. dotted 
line: Eq. (fTf)) . 



Network of A'' patches. Our results can be general- 
ized to a connected network of N patches with migration 
rate /Xy between any two connected patches i and j. We 



liMij, where My 



M. 



3^ 



1. For n = 0,T 



is given by Eq. ([Ij. For /i oo, the rescaled population 
size of each patch becomes Xi — Q/N, where the total 
size of the meta-population Q = Xi slowly varies in 
time. We can perform a canonical transformation from 
xn to the new coordinate Q, keeping xi,X2, ■ ■ ■ ,xn-i 
the same. This procedure [20| brings about a slow one- 
population Hamiltonian HgiowiQ, P), with P = p^^^ and 



effective carrying capacity k from Eq. ([2]) 

1 



1 



2k 



2P 



(18) 



This immediately yields the announced result ([3]). 

How does synchronization of the most probable local 
extinction paths for small /i work for the N patches? 
For /i — > the WKB-instanton is described by Xi{t) = 
Kiq{t — Ti), Pi{t) = p(t — Ti), where Ti are constants, 
i = 1, 2, N, and functions q{t) and p{t) are defined in 
Eq. (fT3|) . This leads to the action 

5(m ^ 0) = 2(1 - In 2) ^ ^ InT^^o/K, (19) 

i 

as announced in Eq. The relative time shifts t;, 

i — 1,2, . . . , N — 1, determining synchronization of the 
local extinction paths, can be found similarly to Eq. (|15|) . 
by minimizing AiS: 



(20) 



AS = — /i max / T)dt, 



where is the subset of indices, corresponding to the 
patches directly connected to the patch i. As in many 
other problems with multi-dimensional instantons j43| . 
the minimization will typically give a unique solution up 
to the overall time shift. Since for the unperturbed local 
instantons Pi{t) is independent of i, we have Xi(^j 0) = 
in Eq. ([20]) . Therefore, AS' < 0, i.e., T is a non-increasing 
function of /i, as in iV = 2 case. Generically, AS* is 
strictly negative, so T decreases with an increase of fi for 
small fi. If all patches have the same carrying capacity, T 
is constant in the WKB regime, up to a pre-exponential 
factor, and corresponds to a single-patch MTE with the 
combined carrying capacity KN IT], |2J|. Finally, the 
validity of the WKB theory demands fin ^ K~^, where 
n is a typical number of connections of a patch (that 
is, a typical node degree of the network), i.e., ~ n, 
whereas the WKB perturbation theory for small fi de- 
mands /xn <C 1. 

In summary, we have developed a quantitative the- 
ory of stochastic extinction of an established meta- 
population where individuals can migrate between dif- 
ferent habitat patches. We have found that, as the 
meta-population goes extinct, local extinction paths be- 
come synchronized already at very small migration rates. 
Finally, we have shown that the MTE of the meta- 
population reaches its maximum for a small but nonzero 
migration rate. B.M. was supported by the Israel Science 
Foundation (Grant No. 408/08), by the US-Israel Bina- 
tional Science Foundation (Grant No. 2008075), and by 
the Michigan Center for Theoretical Physics. 
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Appendix 

We present here derivations of several results outlined in 
the main text. 



Effective carrying capacity of the meta-population in 
the fast-migration limit 

Here we provide more detail on the deterministic dy- 
namics in the fast-migration limit. The rate equations 
for N patches with local carrying capacities Ki — Kni 
are the following: 



'■ /ki + fi^^ {xj - Xi) , (Al) 



for i = 1, . . . , A^. Here li is the set of indices of patches 
connected to patch z, and Mij = Mji. For /i — > oo we 
obtain 

M,, {x, ^ X,) = 0, i = l,...,N. (A2) 



Equations (|A2|) are solved by Xi = Xj = x for all i and j. 
Now, summing up Eqs. (|Aip and replacing Xi by X/N = 
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X, where X = Xi is the total meta-population size, 
we obtain: 



X = X - X^/R, 



(A3) 



Therefore, in the fast- migration hmit, the total meta- 
population size evolves as if the population were oc- 
cupying a single patch with effective carrying capacity 
K = Kk, as announced in Eq. (2) of the main text. 



Slow-migration limit in WKB regime 

General case: N patches 

To calculate the correction to WKB action in the first 
order in /x, described by Eq. (20) of the main text, we 
need to evaluate the following integral: 



I{r) = q{t) 



„p(t-T)-p(i) _ Y 



dt, (A4) 



so that 



AS = — /i max 2, i^iMijl (tj — 



The integral /(r) can be evaluated analytically, using 
the explicit one-patch instanton solution, see Eq. (13) of 
the main text. Upon the change of variables ^ = e*, we 
obtain 



/(r) = 2(e- - 1) 



2 [(e^ - l)ln2 - r] 
-2 ■ 



(A5) 



The first-order correction to the action becomes, there- 
fore, 



AS = u min > 
{n} ^ 



2kjMjj [(e^^'"^' - 1) In 2 - t,- + t,' 

2 - e^3-^' 



(A6) 



The minimization of this sum with respect to the set 
of time differences — tj will generically give a unique 
solution, up to an arbitrary overall time shift tq . 



Example: Two patches 

For iV = 2 we can put ki = 1, K2 = k and M12 = 1. 
Setting, without loss of generality, tq = 0, we obtain from 
Eq. (fXel) : 



AiS = 2fi min 



(e^ - l)ln2-r 



{e-^ - I)ln2 + T 
2 - e-^ 



(A7) 



The minimization of this expression can be easily per- 
formed numerically. For k = 0.25, one obtains AS = 
— 0.5455 ... /i. This asymptotic is shown as the dashed 
line in Fig. 2 of the main text. 



Fast-migration limit in WKB regime 

In the limit of large fi it is convenient to rescale time 
by /i and present the Hamiltonian as iJ = Ha + eHi with 
e = l/^t < 1: 

^"=^E E N i^"''" 1) +-^. (^"'""^ - 1)] ' 



N 



i/i = E-^(^"-i) + ir(^"'"-i) 



It is the "migration Hamiltonian" Hq that dominates the 
dynamics in this limit, whereas the branching and anni- 
hilation terms are small perturbations. To exploit the 
expected time-scale separation, we perform a canonical 
transformation from x and p to qi,Q and Pi^Pq, where 
i = 2,3,...,A^: 

Q = X^^Xj, qi = Xi, 

j 

Pi = Pq, Pi = Pq+ Pi- 

This canonical transformation is motivated by our expec- 
tation that, at small e, the total population size Q and 
the corresponding conjugate momentum Pq are slowly 
varying quantities. 

In the new variables, the Hamiltonian acquires the 
form 

H{{q^},QdPq.},PQ) = ^o(ta},Q,K}) 

+ eHl{{1^},QAPq.},PQlm 



for z = 2, 3, N , where 
^0 = ^ E E ^^^^ (e 



P.-P^ 



1) + q, 



-1)] 



E^i^- 

jell 



N 

E' 



-,-p, 



-1 



(A9) 



N 



Hi - 



Q - * 

2kT 



-2Po 



1) 



(AlO) 
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Notably, the migration Hamiltonian HQ{{qi} , Q , {pq.}) 
does not depend on Pq. This reflects the simple fact 
that, in the absence of branchings and annihilations, the 
total population size Q is conserved. For brevity we shall 
omit tildes in H,Hq and Hi in the following. 

The Hamiltonian equations of motion associated with 
are: 

q^ = ap.ffo(q,g,P)+eap.iIi(q,Q,P,PQ), (All) 
P, = -5g,iIo(q,g,P)-ea,,J?i(q,Q,P,FQ),(Al2) 
Q = e9pi/i(q,g,P,PQ), (A13) 
Pq = -dQHoiq,Q,P) - edQHi{q,Q,F,PQ),iAU) 

where q = (92, 93, gw) and P = {P2, P3, ■-, Pn)- We 
note that no approximations have been made so far, and 
Eqs. (|Alll) - (|A12p are valid for any e. 

The perturbation scheme that we have developed for 
Eqs. (IA11[) - (|A14I) assumes the following time-scale sep- 
aration scenario. At e <C 1 the total population size Q 
varies in time on a long time scale of order e^^ ^ 1. The 
local population sizes qi quickly (on a time scale of order 
unity) adjust to the instantaneous value of Q. Moreover, 
at given Q, the dynamics of the local population sizes are 
essentially deterministic. Having this scenario in mind, 
we can look for the solution as qi — qi{Q) + eq-^^ and 

Pi — eP- ^\ where qi{Q) is the fixed point of Eq. (jAlip 
with e = for an instantaneous value of Q, see Eq. (IA17[) 
below. The smallness of Pi reflects the expectation that 
the fast variables evolve almost deterministically. In this 
case Hq in Eq. (|A9p can be linearized with respect to Pi , 
and we obtain: 



iIo(q,Q,P)=eP'''-h(q,Q)+0(£2), 

hk = ^ Mkj {qj - qk) 

+ MijSkj iQ-qk 
3&I1 \ 




(A15) 
(A16) 

2,3,. ..,A^, 



where Skj is the Kronecker delta. It is seen from 
Eq. (jA15|) that Hq is actually of the first order in e, 
rather than of the zeroth order as could have been naively 
expected from Eq. (jASp . In the leading order in e 
Eqs. (|Alip and (|A12j) become 

= h,{q,Q), (A17) 

= a,,p(i) •h(q,g) + 5,.iii(q,g,o,PQ). (A18) 



Solving the set of algebraic equations (IA17[) for qi, and 
using Eqs. (jA15[) and (IA16p . we obtain 



q,{Q)^Q/N, * = 2,3,...,iV, 



(A19) 



in agreement with our mean-field results in the fast- 
migration limit, see Sec. 1. 

Now let us take the total derivative of Eq. (|A17|) with 
respect to Q: 



dP'-^' ■ h(q, Q) dqj dl"^'' ■ h(q, Q) 



(1) 



dqj dQ 



dQ 



0, (A20) 



and combine this result with Eq. (jA18[) multiplied by 
dqi/dQ and summed over i from 2 to iV. We obtain 



aiIi(q,Q,0,PQ) dqj _ 9PW • h(q,Q) 



dqj 



dQ 



dQ 



(A21) 



where the summation convention is used. Using 
Eqs. (|A15I ) and (|A21|) in Eq. (IA14|) . and restoring the 
original time, we can rewrite Eqs. ()A13p and (|A14p as 

Q = dpHi{q{Q),Q,0,PQ) = dpHsioAQ,PQ), (A22) 

p = -a,^.Hi(q(o),g,o,PQ)^-aQi/i(q(Q),go,PQ), 



— —dQHslowiQ,PQ), 



where we have introduced the slow Hamiltonian 



(A23) 



Hslow {Q, P) 



Hi{qiQ),Q,0,PQ) 

N 



i=2 



N 



Q 

N 



- 1) 



Q^ 



= Q(e^«-l) + ^(e-2^«-l), (A24) 

where R is defined in Eq. (|A3p . This concludes our 
derivation of the effective one-population Hamiltonian 
[see Eq. (18) of the main text] which describes slow large 
fluctuations of the total meta-population size. 



